variables

%s = tf('s');
%testt = (2.1)/(0.1*s^2 + 0.2*s +1); %transfer function

stepsize=0.01;
tend=30; 
ts=(0:stepsize:tend)'; %sampling times
N=length(ts);




%       -1,-2,-3
ycoeff=[0.15;0.27;-0.43];
ucoeff=[0.73;3;7];

yl=length(ycoeff);
ul=length(ucoeff);

coeff=[ycoeff;ucoeff];



omega=5; 
testinp=@(t) sin(omega*t.*t); %input function
u=testinp(ts); %inputs at times




testt=tf(ucoeff',[1;-ycoeff]',stepsize,'Variable','z^-1');
% SYS = TF(NUM, DEN, TS, 'PropertyName1', 'PropertyValue1')
%creates a discrete-time transfer function with
%sample time TS (set TS=-1 if the sample time is undetermined), 
%with numerator(s) NUM and denominator(s) DEN.  
%The output SYS is a TF object.

ym=lsim(testt,u,ts);



y=ym;


%yl+1

firstentry=max(yl+1,ul);

ys=y(firstentry:N);


%[y(yl:N-1),y(yl-1:N-2),y(yl-2:N-3)...]

S = zeros(N-firstentry+1,length(coeff));



for k=1:yl
    S(:,k)=y(firstentry-k:N-k);
end

for k=0:(ul-1)
    S(:,k+1+yl)=u(firstentry-k:N-k);
end




%yp=lsim(testt,u,ts); %y_pure

disturbance=random('norm',0,ones(size(ys))*0.02);

%später noch durch transfer function jagen??

yd=ys+disturbance;

plot(ts(firstentry:N),[ys,yd])




phat=S\yd;


